Electron-electron interactions and the phase diagram of a graphene bilayer 
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We study the effects of long and short-range electron-electron interactions in a graphene bilayer. 
Using a variational wavefunction technique we show that in the presence of long-range Coulomb 
interactions the clean bilayer is always unstable to electron and hole pocket formation with a finite 
ferromagnetic polarization. Furthermore, we argue that short-range electron-electron interactions 
lead to a staggered orientation of the ordered ferromagnetic moment in each layer (that is, c-axis 
antiferromagnetism) . We also comment on the effects of doping and trigonal distortions of the 
electronic bands. 
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I. INTRODUCTION 

The recent developments in the field of carbon physics, 
where a few layers or even single layers of graphene have 
been isolated, have shown that the physics of these sys- 
tems is unconventional from the point of view of tradi- 
tional semiconductor and Fermi- liquid physics 1 ' 2 . The 
electronic dispersion of graphene close to the two K- 
points of the Brillouin zone can be written as 3 -: E±(p) = 
±vp\p\, where vf is Dirac- Fermi velocity (this expres- 
sion is valid for two-dimensional momentum p = (p x ,p y ) 
such that |p| < A where A is a momentum cut-off of 
the order of the inverse of the lattice spacing a). This 
dispersion relation is identical to the one of Dirac elec- 
trons with "speed of light" given by vf ■ In this case the 
electron effective mass, m* , is zero, and the density of 
states vanishes at the K-point. The vanishing of the ef- 
fective mass, the interplay of interactions, disorder, and 
extended defects, lead to anomalous behavior in many 
physical properties-^. 

The capability of experimentally controlling the num- 
ber of graphene layers opens up the field for the study 
of the effect of interlayer coupling in a strongly inter- 
acting two-dimensional system. Interlayer coupling is a 
controversial topic in the graphite literature where the 
precise nature of the coupling between graphene planes is 
unsettled^. Another important issue in carbon research 
has to do with the weak ferromagnetism in highly disor- 
dered graphite that have been observed in experiments 8 
but is still a theoretically open problem*^. 

It is well-known that the low-density electron gas with 
long-range Coulomb interactions in two and three dimen- 
sions is unstable toward a ferromagnetic state. The origi- 
nal argument due to Bloch relies on a variational calcula- 
tion of the ground state energy^. Recently this approach 
was used to look for a possible ferromagnetic instability 
in a single layer of graphene^. The parameter that con- 
trols the relative strength between kinetic and Coulomb 
energies is the dimensionless coupling, g = e 2 e^ 1 /vf , 
(Ti = 1) where e is the electric charge (e 2 = 14.4 eV A), 



and eo is the graphene dielectric constant (eo ~ 1). In 
that case, ferromagnetism is only found for values of g 
larger than a critical value, g c ~ 5.3, which is larger than 
its estimated value in graphene (g ~ 2.1). An analysis 
based on short-range interactions seems to confirm this 
picture . 

In this paper we use a similar variational technique to 
study a clean graphene bilayer where we include the hop- 
ping between graphene planes. Unlike the case of a single 
layer, we find that the bilayer is always unstable toward 
a ferromagnetic state with formation of electron and hole 
pockets with a polarization of the order of 1CU 6 to 1CU° 
electrons per carbon. This result may have direct impli- 
cations for the interpretation of the magneto-transport 
data in graphitic devices^. 

The paper is organized as follows: In Section |TT] the 
model is introduced. In Section lllll we explain the vari- 
ational calculation and present the phase diagram. The 
influence of other hopping parameters on the instability 
are discussed in Section [iVl Section [V] includes the re- 
sults for the low-energy susceptibilities and a discussion 
of short range interactions. The conclusions of the paper 
are to be found in Sect ion IVT1 We also include appendices 
with some mathematical details. 

II. THE MODEL 

The lattice structure for the bilayer which is just one 
unit cell of graphite is depicted in Fig. [TJ For simplic- 
ity we model the system by the nearest neighbor tight- 
binding Hamiltonian: 

Ht.b. = -t (cA^m^Bi^a+h.C.) 

<m,n> 

Z,CT 

m,er 

where c m a (c\. m ff ) annihilates (creates) an electron 
on site m of the sublattice a (a — A, B) of plane i 



2 





FIG. 1: (Color online) Lattice structure of the bilayer. The 
A-sublattices are indicated by the darker spheres. 



FIG. 2: Band dispersions near the K-points in the bilayer. 
Bands are labeled by the numbers 1 — 4 as in the text. 



(i = 1, 2), with spin a (a =f, |), t (t 3 eV) is the in- 
plane hopping energy and t± (t± ~ 0.35 eV in graphite 6 ) 
is the hopping energy between atom Ai and atom A2 (see 
Fig. [T]) . A similar tight- binding Hamiltonian for graphite 
and the single graphene layer was studied long ago by 
Wallace^. At low energies and long wavelengths, the ki- 
netic Hamiltonian can be expanded around the K (K') 
points in the Brillouin zone. The resulting Hamiltonian 
can be written as: 



a 



(2) 



and V 



(p,ai,a, a) 

^p,rA2,cr,a' p.rB2-(7.a/' 



Here, 



with Q denoting 

(J J 

V u p,rAi,u,ai °p,rSi,cr,a' 

c p a t a a ( c p.a,a,a) creates (annihilates) an electron with 
momentum p, on sublattice on (a = A,B) of plane 
i (i = 1,2), with spin a (a =T;i) at the K-point a 
(a = 1,2) in the Brillouin zone, and 



i = 1,2,3,4. The non-interacting ground state has de- 
generacy 4 per momentum, per plane, due to the SU(2) 
spin rotation symmetry and the Z2 real space sublat- 
tice exchange symmetry (at low energies this symme- 
try becomes SU(2) for the continuous rotation of the 
K and K' states in momentum space), and occupation 
at half-filling, given by: rai )ffj0 (p) = 0, n 2 ,a,a{p) = 1, 
n 3,cr.a{p) — and n4 )CT)a (p) = 1. Hence, the presence of 
t± does not mix the spins or the K-points. However, the 
two Dirac cones transform into vertex touching hyper- 
bolae, and for p <C t± the electrons acquire an effective 
mass, m* rs tj_/2. 

The Coulomb interaction in the bilayer is conve- 
niently written in terms of the Fourier components 
of symmetric and anti-symmetric combinations of the 
layer densities, p±(q) = pi(q) ± /0 2 (q), where p;(q) = 
Ek,a,<T,a c k+ q , Ql ,<T,a c k,a I ,^a- The Coulomb term reads: 
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(3) 



is the kinetic energy matrix where 4>(p) = tan -1 (p y /p x )- 
We have set vp = 1 = A, so that the energy is mea- 
sured in units of the in-plane hopping, t , and distance is 
measured in units of carbon-carbon distance a (a « 1.42 

A)ft. 

The kinetic term can be diagonalized by a unitary 
transformation: ^ p = A^p^p, where M p is given in Ap- 



pendix [XJ Then Hkin 
where the four energy bands are given by 



E PJ >,, £ j(p) $ w>,« $ w>.«' 



Ei( P ) 
E 2 ( P ) 
EM 
Ei( P ) 



-t x /2 + E(p), 
t x /2-E(p), 
t x /2 + E(p), 
-t x /2 - E{p) , 



where E(p) = \Jt\ /4 + p 2 . The bands are sketched 
in Fig. [21 Any state of the system can be labeled in 
terms of the occupation of each band, ^i )(7 ,o(p)j with 



Hl= 2S^2Y1 P«(q)V r a (q)p a (-q) 



(4) 



q#0 a=± 



where V±(q) = 2ire 2 (l ± e~ qd )/2eoq, S is the area of the 
system, and d is the interplane distance (d ~ 2.4a rs 
3.35A, Ad « 3.7). We are going to show that in the 
presence of Eq. ([4]) the non-interacting ground state is 
unstable. To perform the calculation it is convenient 
to express the density operators in the diagonal ba- 
sis: p±(q) =E p , !J - j!a $p +q! ,,, i jy(p+q,p)*p, J >,« and 

write the exchange energy associated with Eq. (01 as: 



E, 



1 r dp d 2 P ' 

2 J (2tt) 2 (2tt) 2 



E 



Xij(p\p)x%(P,p')ni,<T,a(p')nj,v, a (p)V a (p' - p). (5) 



The definitions of the matrices an d some more details 
about the exchange interaction for Bloch electrons and 
Eq. §5§ are given in Appendix [Bl 
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III. VARIATIONAL CALCULATION AND 
PHASE DIAGRAM 

Consider the half-filled case with a variational state 
with one electron pocket in the spin up channel and one 
hole pocket in the spin down channel at each K-point: 
ni,T,a(p) = @{Q—p), ni,4. )0 (p) = 0, n 2 ,t, a (p) = 1, and 
n 2,i,a{p) = 1 — ®{Q—p)i where Q, the size of the pocket, 
is a variational parameter (in what follows we assume 
Q <C t± and hence the occupations of bands 3 and 4 
are not affected). Pictures of the non- interacting ground 
state and the trial state are shown in Fig. [3^ and [SJd. 
Notice that the size of the pocket is the same in different 
channels because of the conservation of the number of 
electrons at half-filling. This state breaks the SU(2), but 
not the Z2 symmetry, and is therefore spin polarized (fer- 
romagnetic). There is a similar state that breaks both 
symmetries and has no net magnetization: an electron 
(hole) pocket in the up (down) spin channel in K-point 
1 and a hole (electron) pocket in the up (down) channel 
in K-point 2. We can show that the spin polarized state 
is lower in energy (see below). 



the kinetic term that is of order Q 4 



rn 



There- 




lb) 



(d) 



z_ _^ z_ _^ /_ _\ z_ _\ 

FIG. 3: (Color online) Sketch of the trial states: (a) 
Half-filled non-interacting ground state, (b) Trial state 
with particle-hole pockets built upon (a), (c) Doped non- 
interacting ground state, (d) Trial state with particle-hole 
pockets built upon (c). 



The change in the kinetic energy per unit area due to 
an electron (or hole) pocket of size Q is given by: 



AE kin _ 1 r (Q 2 +tj/4) 3 / 2 tj 

S 2vr L 3 24 



t±Q 2 



.91 

&7rt 1 



( 6 ) 

up to order Q 4 . The expressions for the change in the ex- 
change energy are cumbersome and details are provided 
in Appendix [Cl We find from Eq. |(SJ), up to the same 
order: 
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Notice that the leading order term in the exchange in- 



teraction is AE ox /S 



-5Q7277T 2 



i 3 / 2 , where 



fore, we have proved that the bilayer is always unstable 
to the formation of polarized electron and hole pockets. 
In contrast to the single graphene plane casein, the to- 
tal energy is negative for small Q. This is due to the 
fact that the exchange with the filled bands is less im- 
portant in this case. In order to calculate the equilib- 
rium size of the pockets we minimize the total energy, 
AE tot (Q) = AE kin (Q) + E ex (Q), with respect to Q and 
find Qmin, that is, the size of the pocket for which the en- 
ergy is minimized. For the parameters in graphene (see 
below) we find that Q m i n ~ 0.05<j^ (<C t±), justifying the 
above expansion. 

Consider the case where the system is initially doped 
with pockets of size Qo ( n — Qo/2). We look for an 
instability by varying the density of electrons and holes 
subject to the constraint of particle conservation. Note 
that the instability can produce one type of carrier (either 
electron or hole) if Qo > Q m i n or two type of carriers 
(electrons and holes) if Qo < Qmin- We can parameterize 
the state with one type of carrier by taking Q 2 = Qq (2 — 
x) and Q 2 = Qq x with < x < 1. For the state with 
two types of carriers we take instead Q 2 = 2Qq + \x\ 
and Q 2 = \x\, with x < 0. The doped non-interacting 
ground state and the trial state with particle-hole pockets 
(x < 0) are shown pictorially in Fig. [3J: and [3Ji. The 
calculation proceeds as before and we find: 

AE 1 

— w —[AE tot (Qt) + A£ tot (Q x ) - 2A£ tot (Q )] 

>(Qo,x). (8) 



-AE,. 



The extra term, A_E cx t r a(Qo: x )i comes from terms that 
cancel out in the undoped case. To leading order in Q 
these terms are given by <?<2o(l ~ s) 2 /(647r) for x > 
and gQoiQ 2 , + 2|a;|)/(647r) when x < 0. In our units we 
have tj_ -C 1 so that, to a first approximation, this contri- 
bution is much smaller than the quartic term in Eq. {7}, 
and it can be neglected. This leaves us with the first line 
in Eq. {S} involving Ai?tot only. The dependence on x is 
implicit through Qf and Q^. Then Eq. (J5J) has the form: 
AE(Q) = -A\Q\ 3 + B\Q\ 4 . Rescaling the Q variable so 
that the minima of the energy in the paramagnetic states 
sits at \Q\ — 1, we have: 



A£(Q) = -|Q| 3 /3 + |Q| 4 /4. 



(9) 



m is the magnetization, which is always dominant over 



Using the scaled variables we see that the system is un- 
stable to small deviations in x from 1 if Qo < 1/2. The 
ferromagnetic state has lower energy than the param- 
agnetic states if Qo 0.7 and the resulting state has 
electron and hole pockets. As a consequence of the first 
order nature of the transition, the system exhibits phase 
coexistence (that can be obtained from a Maxwell con- 
struction, not shown in the Fig. 2]), and hysteresis in 
physical properties such as magneto-transport, around 
the critical line. In this region, the system shows a ten- 
dency towards electronic phase separation^, frustrated 
by electrostatic effects. As the charge densities involved 
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are rather low (see below) we cannot exclude the forma- 
tion of large domains of the different phases. The phase 
diagram for t± = 0.05 is shown in Fig. [4] as well as a plot 
of AE(Qq,x) for some typical cases of Qo- 
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FIG. 4: (Color online) Left: Phase diagram of the graphene 
bilayer as a function of electron density away from half- 
filling, n (electrons per carbon), and coupling strength, g — 
e 2 /(eoVF), with t± = 0.05. Inset: AE as a function of x (as 
defined in the text) in the paramagnetic (A), critical (B), and 
ferromagnetic (C) regions of the phase diagram. 

In the previous calculation we have not included the 
exchange interaction between different K-points in the 
Brillouin zone. In that case the spin polarized state 
that breaks SU(2) is degenerate with the state that 
breaks both SU(2) and Z2. The difference between 
the states is how the pockets are assigned to the spins 
and the K-points. By including exchange between K- 
points in Eq. ([5]) we find that there is small energy 
difference between the states favoring a state with a 
net ferromagnetism but which retains the Z2 symme- 
try. Quite generally this is the case since the elements of 
X?j(p' >p)X%(p>p') are all positive. A direct calculation 
using Eq. and taking p and p' to lie at nearest neigh- 
boring K-points confirms this picture. One then finds a 
very small energy difference of order oc Q 4 , and hence 
other effects can be important in determining the actual 
ground state. There is also a correction to the Q 4 /t±- 
term in Eq. |[7J) that changes the position of the optimal 
value of Q m in by a small amount. 

In order to compare with experiments it is interest- 
ing to estimate the total magnetization in the polarized 
state for the case of the undoped graphene bilayer. We 
estimate the cut-off using a Debye approximation so the 
the number of states is conserved in the Brillouin zone: 
A 2 = 2tt/A u = 47r/(\/27a 2 ), where A u is the area of the 
real space unit cell. Restoring the units, we set vf = 
37 a/2 « 10 6 m/s, and t± = ji/(v F A) ~ 5.2 x 10~ 2 , 
where 71 » 0.37 eV is the typical graphite value 6 . Hence, 
for two pockets of size Q the density of electrons per car- 
bon is approximately n = Q 2 /4 ps 1.6 x 10~ 6 (tj_ ~ 0.05 
and Q ~ 0.05t±), and therefore, the magnetization per 
carbon is m m 10~ 6 — 10~ 5 /is (/jb is the Bohr magneton) . 
These number are, of course, very approximate because 
the value of the microscopic parameters do not need to be 



the same as in graphite and the presence of a cut-off intro- 
duces further uncertainty. In any event, the magnetized 
state of the graphene bilayer shows very weak ferromag- 
netism. A direct experimental consequence of our calcu- 
lation is that the bilayer has two species of electrons (elec- 
trons and holes) and therefore they should contribute to 
the Hall resistivity at small magnetic fields, B. In par- 
ticular, it is easy to show that the magneto-resistance at 
small magnetic fields acquires a B 2 dependence 1 ^. 



IV. OTHER HOPPING PARAMETERS 

We would like to comment on other effects that we have 
not considered in the previous calculation. In terms of 
the Slonczewski-Weiss-McClure model for graphit o 15 i 16 
our model includes only the parameters 70 (t) and 71 
(ijj but not 73 and 74. On the one hand, 74 introduces 
an electron-hole asymmetry by changing the curvature 
of the bands, but the bands remain parabolic near half 
filling. On the other hand, 73 introduces a trigonal dis- 
tortion which restores a linear dispersion at low energies. 
To estimate the effects of 73 we use the effective low- 
energy model that can be derived from the extension of 
Eq. §3§ to include 73 and projecting onto the two bands 
that are closest to the Fermi surface (see Ref. [I?} and 
Appendix [El for details). The effective kinetic energy ma- 
trix is then: 



/C(p) 



p2 / 

t± I 



V3P 



e l<t 
e-** 



(10) 



where V3 — 73/70, with energy eigenvalues given by: 

E±(p) = ±pVp 2 + M±) 2 + 2pv 3 t± cos[30(p)]/*j_.(ll) 

It is easy to see that the crossover from linear to quadratic 
dispersion takes place at a momentum p C ross ~ ^3 1 j_ . 
If Across *C Qmin (i>3 <C 0.1) the previous calculation is 
still valid since the dispersion remains parabolic at the 
scale of the instability. Nevertheless, if one uses the val- 
ues of the parameters in the graphite literature 6 -, namely, 
v 3 = 73/70 ~ 0.1 we conclude that j? C ross ~ Qmin and the 
trigonal distortion may become important. We should 
remark, once again, that it is not guaranteed that the 
value of the parameters in graphite are the same as in 
the bilayer (hence, one should take the numbers here 
with a grain of salt). At small energies, the spectrum 
in Eq. (fTU|) can be described in terms of one Dirac cone 
at p — 0, and three asymmetric cones at p CO nc = v 3 t± 
(the direction of the cones are such that cos(3</>) = — 1). 
This situation maps onto the single graphene plane case 10 
with ellipsoidal pockets instead of circular, and a small 
renormalized Fermi velocity, v* F m v^vf (<C vf)- There- 
fore the dimensionless coupling strength is g* = g/i)s s» 
20, which is much larger than the critical coupling for 
ferromagnetism 10 . Then, in this case the transition from 
paramagnetism is into a fully polarized state with effec- 
tive bandwidth of order of v^t^, leading to a polarization 
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of the order of 10~ 6 electrons per carbon, which is of the 
same order of magnitude of polarization found without 
73 . Unfortunately this argument is not rigorous since the 
exchange with the filled bands also contributes, requiring 
a more detailed study. 



V. SHORT RANGE INTERACTIONS 

We consider the effects of short-range electron-electron 
interactions which, for simplicity, we describe by an on- 
site Hubbard interaction, U. It turns out that this inter- 
action favors an antiferromagnetic ordering. In order to 
quantify the tendency towards this phase, we calculate 
the associated susceptibility and present a simple mean 
field argument. 



A. Electronic susceptibility. 

Using the basis defined in Eq. (fTU|) with V3 — (the 
procedure outlined in Appendix [E| the wavefunctions 
corresponding to the two bands closest to the Fermi level 
can be written approximately as: 



*±(k) 



I / e -i<Hk) 
71 V ±e 40(k) 



(12) 



so that the states are mostly localized at the sites with- 
out a corresponding atom in the neighboring layer. For- 
mally, this two component wavefunction is equivalent to 



the spinor defined in the analysis of a single graphene 
plane except that the angle is doubled. Restricting the 
calculation to this subspace, we can write the bare sus- 
ceptibility as a 2 x 2 matrix (Ref. HT 



xo(q,w) 



XD(q,w) XND(q,^) 
XND(q,w) XD(q,w) 



(13) 



Here xp denotes the response in the same plane as the 
source and xnd response in the opposite plane. The ran- 
dom phase approximation (RPA) susceptibility, assum- 
ing an on-site Hubbard interaction U, is: 



XR PA (q,^) = xo(q,w)[i - c/xo(q,w)]" 



(14) 



This expression becomes simpler when decomposed into 
a contribution symmetric in the two sublattices and an- 
other antisymmetric: 



XFM(q,^) = xp(q^) + XNp(q,^) 

4 1 - U[xn(q,u) + XNp(q,w)] 

XAFM(q,^) = xp(q^) - XNp(q,^) 

4 1 - C/[xp(q,w) -XNp(q,w)] 



(15) 



The symmetric susceptibility gives the response of the 
system to a magnetic field which is the same in the two 
sublattices (note that we are neglecting the influence of 
the sites where the states have zero weight of k = 0), 
and induces a ferromagnetic ordering. The antisymmet- 
ric response leads to antiferromagnetic ordering. The 
susceptibilities can be written as: 



XD(q,w) 

XNp(q,w) 



tx 



tx 
4tt 



^ ,\ dk r d ± 

4vr J J 2tt 

kdk / 
i) Jo 



1 



1 



d<t> 

2^ 



wt ± -|k + q/2| 2 -|k- 

COS (20 k+q /2,k-q/2) 



q/2| 2 ^ ± + |k + q/2| 2 + |k-q/2| 2 

COS (20 k+q /2 :k _ q / 2 ) 



ujt x - |k + q/2| 2 - |k - q/2| 2 ut x + |k + q/2| 2 + |k - q/2| 5 



(16) 



where </>k+q/2,k-q/2 is the angle between k + q/2 and k — q/2, and: 



COS 2 (<^k+q/2,k-q/2) 



(k + q/2).(k-q/2) 
|k + q/2||k-q/2| 



(k 2 - q 2 /4) 2 



{k 2 + q 2 /4) 2 -k 2 q 2 cos 2 (0) 



(17) 



The only dependence on the angle 8 between the vectors k and q of the expressions in Eq. (fl6| is (after using the 
double angle formula) through the cosine in Eq. (fTT|) . Averaging over angles, we obtain: 



(COS 2 (0k+q/2,k-q/2)) = 



|k| 2 -|q| 2 /4 



|q| 2 /4 



(18) 



Inserting this expression into Eq. ([To]) it is a simple task to perform the remaining one-dimensional integral. Intro- 
ducing Xfm = Xp (q ; w ) + Xnp (q, and Xafm = Xp (Qj w ) — Xnp (q, w) we can extract the leading dependence on the 
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cut-off A of the susceptibilities, and we finally obtain: 



R- e XAFM(q^) 



-Ef 
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s=±l 
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s=±l 
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RexFM(q^) > 



Im XAF M (q^) = y 6, (k|-^-jsign(w)-Imx FM (q^)- 



sign(w) , 



(19) 



Hence, setting q — 0, the antiferromagnetic susceptibility 
diverges logarithmically with the cut-off for any finite 
frequency to. The logarithmic dependence implies the 
existence of an instability for any positive value of the 
interaction U. Alternatively, we can show the existence 
of this instability by a direct calculation of the correlation 
energy gained by polarizing the system. 

It is worth noting that the divergence obtained here, 
and the related marginal behavior of a local interaction 
can be obtained from the same power counting arguments 
used for the analysis of two dimensional interacting elec- 
trons near a van Hove singularity 19 . If we include next 
nearest neighbor couplings through the parameter 73 , as 
discussed in Section [IV] the low energy bands can be 
described by an effective Dirac equation. The screen- 
ing of the long range Coulomb interaction vanishes, and 
the corresponding susceptibility can also be calculated 
analytically^!. 

Note also that the polarization function is simply re- 
lated to the susceptibilities^ 2 - above by n° = — x°. This 
allows one to get the screening properties within the 
RPA easily. In particular, for the mode that is sym- 
metric in the layer densities which originally had the 
long- wavelength l/|q| singularity the static (uj = 0) 
RPA screening cuts off the singularity by taking l/|q| — > 
l/(|q| + <7tf)- From Eq. JTS]) we find that the Thomas- 
Fermi screening wavevector <7tf oc t±. This is in 
agreement with what one expects from the usual two- 
dimensional electron gas where the screening wave- vector 
is oc m* independently of the density of camer a 23 ! 24 . 



B. Mean- field approach 



where the ± label different spin orientations, 
bands are then: 



The four 



e(k) - ±1 



' 2fc 2 + t\ + 2A 2 ± t ± y/4 k 2 + i\ 



(21) 



The opening of the gap lowers the kinetic energy of the 
system. We estimate this by performing the integral up 
to k = A = 1, then for t± = .05 



SK 



kin 



-2 [1.9043- y ln(A)]A 2 + 0(A 3 ), (22) 



is the change in the kinetic energy per unit cell due to 
A. The connection between the average magnetization 
M = I < rij^ > — < rij^i > I and the mean field is A = 
UM/2, where U is the strength of the on-site Hubbard 
interaction. The energy price one must pay per unit cell 
for having doubly occupied sites is 



SE, 



UM 1 



U 



;A 2 



(23) 



Because of the logarithm in Eq. (f2"2")l a small anti- 
ferromagnetic distortion is always favorable. Assuming 
that A is small, the mean-field solution is: 



»MF 



-PI"" 



2/U- 1.9043 \ 
tZf2 J' 



(24) 



and hence A is exponentially suppressed unless U is of the 
order of A. Other variations of the mean-field in Eq. (|20|) 
give similar results. Thus within the mean-field approxi- 
mation, the anti-ferromagnetism found here is very weak 
unless the interaction is strong. 



Alternatively we can explain the diverging susceptibil- 
ity with a simple mean-field approach. We introduce a 
staggered mean-field A into the Hamiltonian according 
to 



iC(p) 



pe 



(p) 



pe 



t± 









tx 




pe 






-i0(p) 



(20) 



pe 1 



±A J 



VI. CONCLUSIONS 

In summary, we have shown that long-range Coulomb 
interactions in a clean graphene bilayer lead to a ground 
state that is magnetically polarized with electron and 
hole pockets. We have determined the phase diagram of 
this model as a function of the coupling strength and dop- 
ing, with a first order phase transition line between the 
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paramagnetic and ferromagnetic states. Around the crit- 
ical line one expects hysteresis effects associated with the 
presence of phase coexistence and/or magnetic domains. 
We have also shown that on-site electron-electron inter- 
actions produce a staggering of the ferromagnetic order 
in the two planes and hence, c-axis antiferromagnetism. 
The introduction of other terms in the Hamiltonian, such 
as trigonal distortions, makes the phase diagram even 
richer, due to the creation of new energy scales. It is 
clear from our studies that graphene bilayers present an 
electronic behavior that is rather different from ordinary 
metals. The study of these systems becomes even more 
relevant given the recent developments in the fabrication 
and control of graphene multi-layers, and their possible 
application in nano-electronics. 
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APPENDIX A: UNITARY TRANSFORMATION 



APPENDIX B: EXCHANGE INTEGRAL FOR 
BLOCH ELECTRONS 

Quite generally the Hartree-Fock energy consists of 
three terms. A kinetic term, a direct charging term, and 
an exchange term. The exchange term is^ 

£ cx = < m, n|U|n, m >, (BI) 

where V is the interaction potential and the sum is over 
occupied states. Expanding in Bloch states V'k.a we get 

E C x = — ^ Hk,a,o-^k',o!',CT'/k ) 'k' ) (B2) 
k,a,cr 
k' ,a' 

where for the unscreened Coulomb interaction in 2D 

^k, a W^, a -W^, tt -(x')4. a (x') 

^ — ■ (B3) 

Let us consider one plane only, then we can write the 
Bloch states in the tight-binding approximation as 



The unitary transformation Ai p that one needs to di- 
agonalize the Hamiltonian in Eq. ([3]) can be written as 
M p = A4i(p)A4 2 A4 3 (p), where 



Mi(p) 



(I 

e - ^( p ) 




\ 



o 







1 

o e^wy 



(Al) 



is a gauge transformation, 



Mi 




(A2) 



forms symmetric/ antisymmetric bands, and 

/ cos[tp(p)} sin[<p(p)] 

KA (tA- ~ sin t^(P)] coe #(p)] 
■M 3 [p)-\ Q Q cos[^(p)] -sin[^(p)] 

\ sin[vj(p)] cos[(/3(p)j 

(A3) 

takes care of the final diagonalizing rotation. Choosing 
tan[2yj(p)] = 2p/t± the rotated Hamiltonian /Cdi ag (p) = 
M p IC(p)M p becomes diagonal: 

/C diag (p) = diag[-t ± /2 - E(p), -f x /2 + E(p), 

t ± /2 + E(p),tj2-E(p)]. (A4) 

Except for the labeling of the states (1 — 4), which is 
just a permutation, this is the unitary transformation we 
need to diagonalize the non-interacting problem. 



4,«W = e lk - TA a k ,aw(x - T A ) 

Ta 

+ ^e lk - TB 6 kiQW (x-T B ), (B4) 

T B 

where w(x) is the localized basis function, Ta (Tb) are 
the lattice vectors of lattice A (B) and a k ,a and 6k, a are 
the functions that generate the Bloch state in question. 
Neglecting the overlap of wave- functions on the A and B 
sites we get approximately 

^(x)Vk',a'(x)~ e i(k '- k) -7-S 

x [ a k. Q a k '. tt ^A(x, k - k') + & k , f A'.a'<Mx, k - k')] ■ 

(B5) 

The functions 4>a and </>b are periodic in the real-space 
lattice and can hence be expanded in components of har- 
monics of the reciprocal lattice {K}. Keeping only the 
leading constant term coming from the K = terms we 
get <j) = 1, in which case 



^k,'k' — [ a k,a a k',a' + ^k,a^k',a'] 



2^e 2 



x [a k , Q /a k ct + &k',a'frk,oJ _ jf/j 



(B6) 



The corrections to this K = term are down by at least a 
factor of |k— k'| 2 as is explained in Ref. [25|. It is possible 
to include higher harmonics in the reciprocal lattice, but 
then the important divergence near k w k' is cut-off by 
K. Moreover, including K ^ we should also include 
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short-range (high energy) physics that is not described 
by the continuum model used here. 

Note also that the expression in Eq. (|B6j) is just what 
one get from a simple Fourier transform if one also in- 
cludes the spinor structure due to the two sub-lattices. 
We apply this to the bilayer where the Coulomb interac- 
tion can be written 

Hl = \J d 2 xd 2 y{^ D (x-y)[pi(x) Pl (y)+p 2 (x)p 2 (y)] 

+ U ND (x - y) [p x (x)p 2 (y) + p x (x)p 2 (y)] } . (B7) 

We Fourier transform this and introduce the symmetric 
and antisymmetric combinations 



P±(q) =pi(q)±p 2 (q). 
to write the interaction in a diagonal form 



(B8) 



The AT's are sums of certain components of the matrices 
^■(p'.p) = Xy(p',P)x?i(P,P')- Note that the elements 
of £ are all greater or equal to zero. 

1. Half-filling 



At half-filling the trial state is characterized by the 
single variational parameter Q. It is straightforward, al- 
beit tedious, to perform the matrix multiplications and 
extract the K's. The results for an electron and a hole 
pocket of size Q is: 



E(p) E(p> 



■ cos 



(0-0'), 



(C2) 



(B9) 



q a— ± 



where V± — ^^-(1 ± e qd )/2. The prime on the q-sum 
denotes that the q = term should be excluded since it is 
canceled by the positive (jcllium) background. In terms 
of the operators that diagonalizes the kinetic terms the 
density operators can be written as 



p+qX ± (p- 



q, p)$ 



(BIO) 



2E(p>) 



oil 



1 - 



(C3) 



K v = 1_P V_ 

2 2 E(p) E(p>) 



(C4) 



where 



x ± (p + q, P) = M 



p+q 



'10 
10 
±1 

vO ±1/ 



M 



(Bll) 



and A^ p is given in Appendix [AJ Using this one can 
easily generalize Eq. (|B2|) and (1B6|) to arrive at Eq. ([5]). 



2 E(p) E(p>) 



1 - 



1 - 



E(p) il E{p' 



[1 



^][l + ^]cos2(^')}. (C5) 



APPENDIX C: EXCHANGE INTEGRAL 

Using U D (q) = 2itg/q and U ND (q) = 2wge-i d /q one 
quite generally get that the change in the exchange en- 
ergy can be written as 



AHi 

< —=- >= 



6(Q - p)0(Q - p') [K»V D + K™V ND ] } . (CI) 



In this appendix we temporarily relabel t±/2 — > t to 
avoid an excessive amount of 2's in the equations. We 
can write 



^ AHi _ gv F A 3 



S 



ND 



(C6) 



where we are measuring all parameters (i.e. Q and t) in 
units of the cut-off A. The C's are given by the integrals 
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L 2 



C ND = 



(I 



Cf = 2Q 3 / / .//(/// / do 

^yND 



2/ cos((/)) 



/o JO 

-Q 2 / xcfe / ydy 
Jo Jo -in 

13 r 1 i- 1 



V< 2 + (Qz) 2 v 7 ^ 1 ^ ly " 

j. „-dA|y-Qx| . 



[1 



xdx 

2 Jo Jo Jo 



ydy ( d^Q 2 \ 



3(20)}, 



■[1 



^WTjQyJ ^/W+jQxY J ly - x l 

)3 /-l /-l -dAQ|y-x| / 

-4" / *<fc / 2^2/ / # - ... (2Q 2 [ 



cos(0) 



4 
[1- 



JO 



[1 



y 



^ + {Qxf il y/t* + {Qyf 



[1 



COs(r/>) 



[1 



cos(2</>) 



(C7) 



From this one can extract the leading and sub-leading 
terms in an expansion in powers of Q, the result is given 
in Eq. ([7]). Some useful expressions for performing the 
expansion are provided in Appendix [Dl 



2. Doped case 

The calculation for the doped system proceeds exactly 
as in the previous case but we must allow for the elec- 
tron and hole pockets to have different size. For each 
electron (or hole) pocket of size Q e (Qh) there is a con- 
tribution like that in Eq. (jClj) . The if 's are half of those 
in Eq. ([05]) . (jC4)) and ([U5]). But Kf is different for holes 
and electrons: 



^le/h 



4 



pp 



2E(p)E(p>) 



cos(0). 



(C8) 



The contributions coming from the ±1/2 are easy to ob- 
tain and can be encoded in a new contribution C ncw in 
Eq. jC6]), where 



Cnew — —Q 2 Ro{Q), 

and i?o is given in Appendix [Dl 



(C9) 



Then one can show that up to 0{Q m ) 
V ^ 5Q 6 



Ro{Q) 
Ri(Q) 



i 



7T 

-( 

G 



8 

'HQ) 



2 

2tt 
~9~ 



3Q 3 
80 



91 

64 

7T 

3 

15 Q 5 
1792 
3Q 2 
16 



35 Q s 



1024 16384 
[ln(2) - ^} Q 

175 g 7 2205 g 9 



55296 1441792 

sg 4 2i g 6 9g 8 

"288" + 4096 + 4096 



Moreover i? (l) = 4/3, R a (l) = 2(2C-l)/3and R 2 (l) 
4(3?r - 7)/27. C » 0.91596 is the Catalan constant. 



APPENDIX E: APPROXIMATE TWO-BAND 
MODELS 

There are two reasons for constructing approximate 
two-band models. Firstly, on physical grounds the high- 
energy bands should not be very important for the low- 
energy properties of the system. Secondly, it is much 
easier to work with 2x2 matrices instead of 4 x 4 matrices. 
In this appendix we derive the low-energy effective model 
by doing degenerate second order perturbation theory. 
The quality of the expansion is good as long as vfP <C t± . 



1. Simple low-energy model 



APPENDIX D: COULOMB INTEGRALS 



Let us define 



Rn(Q) 



dx I dy I d<p 



xy cos(n0) 



\jQ 2 x 2 + y 2 — 2Qxy cos(0) 



(Dl) 



If we transform the Hamiltonian matrix in Eq. ([3]) by 
taking the symmetric and anti-symmetric combinations 
of the first and third rows (and columns) we can write 
Wkin — Ho +Hi, where 



Ti.o 



'-t_L 0\ 



ii 

0/ 



(El) 



10 





( 


pe i<f>(p) 





p e -i<t>{p) \ 


1 


pe-iHp) 





pe -i<t>{p) 















_p e -#(p) 




^ pe 1 ^ 





-p e -i<i>(p) 


o ) 



(E2) 

Performing second order perturbation theory one finds 
an effective Hamiltonian matrix, Tt e g = —n[(l/'Ho)Ti.i 1 
where: 



Hcff 



El 
t± 



/0 \ 

e~ 2l<t >^ 



\q e 2#(p) 0/ 



(E3) 



The low-energy spinors are then given by: 



\ 

e -i<t>(p) 





(E4) 



and the corresponding energies are ±p 2 /t±. If we add 
the contribution from 73 (which is already diagonal in 
this basis) and only keep the B-atom components we im- 
mediately arrive at Eq. (p~0|) . 

In fact it is easier to see the existence of the ex- 
change instability in this basis. Working in this subspace 
we again get an expression like that in Eq. (|C1[) with 
Kf = 0, Kf D = - cos(20), K§ = 1 and K^ D = cos(20). 
If one further neglects the difference between V D and 
1/ ND the resulting change in the exchange energy can be 
expressed with the help of the functions defined in Ap- 
pendix [D] Explicitly the leading term in the exchange 
energy is cx -Q 3 [R (l) + R 2 (l)] +2Q 2 R 2 (Q) ~ -8Q 3 /27 
in agreement with the result in Eq. ([7]). 



2. More general low-energy model 

A more general Hamiltonian model for the low energy 
physics is given by^: 



Ho = 



( 




t± 


Vipe~^\ 







Vipe-^ 


v 3 pe^ 


t± 


Vipe 1 ^ 





Vppe'^ 


\ v^pe^ 




v F pe irt > 


/ 



(E5) 



One can perform a unitary transformation so that H.2 = 
M^HoM, where M = A^i(p)7W 2 and M\{-p) and M2 
are given in Eq. (lAlj) and (jA2|) . We may then separate 
the transformed Hamiltonian into three parts, H.2 = /Co + 
fCi + IC2, with: 



/Co 



IC2 = 



ft ± 0\ 

-t x 

\o 0/ 



v 3 pcos(3cf>) 



iv 3 psm(3cj)) 

( (v F + v 4 )p 
(v F + Vi)p 










(E6) 



(E7) 





< ° 

[Vp — V4)p r 



With this decomposition it is easy to find the approxi- 
mate eigenstates and eigenvalues for vpp <C t± . 



For the high-energy states one can use the simple non-degenerate perturbation theory. The eigenvalues are given 
by E 3 = t± + (vf + V4) 2 p 2 /t± and E4 = —t± — (vp — V4) 2 p 2 /t±. It is also straightforward to obtain the corresponding 
states. For the low-energy sector the second order perturbation result (from two IC2 and one /Co) can give a term 
which is of the same order as that of fCi. Thus, we must use degenerate perturbation theory. The usual manipulations 
then given the Hamiltonian matrix in the low energy subspace as /Ci ow = fC\ — ^Viil/ JCo)ViJC 2 , where V\ is the 
projection out of the low-energy subspace, explicitly 



/Clow — 



/o 

«3pcos(3(/>) 
' 
yO iv 3 psm(Z§) 

and the corresponding eigenvalues are 





-W3psin(30) 


-v 3 p cos(30), 



2vpV4p 2 

t± 



± \ (v 3P ) 2 + 



(vp + vpP 2 
t± 



p_ 

tx 



/0 

-(«f + v 4 


Vo 





2 




- 2(«aP) 



(vp + v\)p 2 



(E9) 



(vp — v±y 



cos(30). 



(E10) 
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